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ABSTRACT 


V 

This  report  describes  the  computer  model  NRFL02  which  has 
been  developed  to  calculate  the  near-field  stratified  turbulent  flow 
driven  by  the  intakes  and  outflows  of  an  ocean  thermal  power  plant. 

A two-dimensional  geometry  is  assumed,  with  the  power  plant  intakes 
and  Jet  outflows  modelled  by  boundary  conditions  on  the  left  boundary 
of  a rectangular  domain.  Horizontal  flow  through  the  right  boundary, 
to  or  from  the  far- field,  is  allowed  by  assuming  a pressure  in  hydro- 
static equilibrium  with  the  ambient  density  distribution.  The  code 
uses  a first-order  closure  model  to  treat  the  stratified  turbulence. 
Advanced  numerical  methods  enable  convergent  and  accurate  solutions 
to  be  obtained  rapidly  and  economically. 

A simple  laboratory  simulation  of  two-dimensional  stratified 
turbulence  driven  by  intakes  and  jet  outflows  is  proposed  in  order  to 

address  the  general  question  of  near-flow  recirculation.)  Such  an 

- — • 

experiment  would  also  provide  vital  data  for  verifying  the  code  NRFL02 
and  its  successors , and  for  tuning  the  turbulence  parameters  to  give 
the  best  agreement.^*  Numerical  results  from  NRFL02  are  presented  for 
the  proposed  experimental  simulation.  They  confirm  our  analysis 
predicting  a stability  limit  on  the  far-field  horizontal  flow,  and 
demonstrate  the  occurrence  of  recirculation  for  fast  flow  cases. 
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1. 


INTRODUCTION 


Estimates  of  the  resource  availability  and  environmental 
impact  of  an  ocean  thermal  power  plant  (OTPP)  will  require  carefully 
designed  computer  models  of  the  flows  driven  in  the  ocean  waters 
immediately  surrounding  the  plant.  This  report  describes  the  computer 
model  NRFL02  developed  to  study  the  stratified  turbulent  motions 
resulting  from  the  intakes  and  outflows  of  an  OTPP.  The  application 
of  NRFL02  to  particular  OTPP  engineering  design  configurations,  such 
as  discussed  by  Trimble  (1975)  and  Douglass  (1975),  will  be  presented 
in  subsequent  reports. 

An  ocean  thermal  power  plant  operating  in  tropical  waters 
uses  warm  surface  water  to  boil  a working  fluid  such  as  ammonia.  The 
high  pressure  gas  generates  power  by  passing  through  a turbine,  and  is 
condensed  at  a low  pressure  using  cold  water  from  the  deep  ocean.  A 
typical  200  MW  power  plant  would  take  in  5 x 104  cu  ft/sec  of  water 
from  near  the  ocean  surface,  at  about  80°F,  and  eject  it  at  about  77°F. 

A similar  volume  flux  would  be  removed  from  the  deep  ocean  (at  depth 
1500  to  U000  ft)  at  about  4o°F,  and  ejected  at  about  U3°F  after  being 
used  to  condense  the  working  fluid.  Most  designers  envisage  floating 
plants,  though  the  early  prototypes  may  be  on  land. 

Our  study  of  the  turbulent  external  flow  near  an  OTPP  has  two 
objectives.  First,  we  wish  to  determine  the  average  inflow  temperatures, 
since  their  difference  constitutes  the  thermal  resource.  For  certain 
designs  and  ocean  environments,  there  will  be  a substantial  reduction  in 
the  average  temperature  of  the  warm  inflow  water,  due  to  turbulent 
mixing  and  to  recirculation  of  some  of  the  outflow  water  from  the  plant. 
Our  second  objective  is  to  obtain  results  from  these  near-field  studies 
which  can  be  used  in  calculating  the  far-field  environmental  impact  of 
OTPP  operation.  A detailed  discussion  of  the  far-field  effects  and  their 
background  in  geophysical  fluid  dynamics  is  given  by  Piacsek,  Toomre,  and 
Roberts  (1975).  Our  calculations  to  be  discussed  here  are  confined  to  the 
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region  near  the  plant,  where  the  flows  can  deviate  substantially  from 
the  horizontal,  and  the  turbulence  is  very  much  stronger  than  in  the 
ambient  ocean.  Far  from  the  OTPP,  the  plant-induced  addition  to  the 
ambient  ocean  currents  is  practically  horizontal,  and  the  plant- 
induced  turbulence  has  decayed  to  negligible  levels. 

In  this  report,  we  describe  a two-dimensional  computer  model 
NRFL02,  which  we  have  used  to  obtain  numerical  results  on  ■the  near-field 
flow.  The  differential  equations  are  presented  in  Section  2,  together 
with  the  first-order  closure  model  used  to  represent  the  stratified 
turbulence.  In  our  formulation,  an  extra  equation  is  solved  for  the 
mean  kinetic  energy  of  the  turbulent  fluctuations.  The  single  turbulent 
diffusivity  is  obtained  from  this  in  a novel  way  which  serves  to 
describe  the  physics  of  stratified  turbulent  transport  and  which  does 
not  allow  the  turbulence  to  decay  to  zero  in  a finite  time.  A further 
innovation  aids  in  obtaining  steady  mean  flows  for  strong  stratification. 

Our  representation  of  the  multiple  inflows  and  outflows 
characteristic  of  an  ocean  thermal  power  plant  is  described  in  Section  3. 
The  .computational  domain  is  rectangular,  with  the  upper  boundary  repre- 
senting the  ocean  surface.  The  left-hand  boundary  is  occupied  by  the 
inlet  and  outlet  ports  of  the  OTPP.  The  finite  computational  domain  is 
intended  to  represent  a semi-finite  portion  of  the  ocean,  and  thus  our 
other  boundaries  (lower  and  right-hand)  must  be  suitably  permeable  and 
with  an  appropriate  thermal  stratification.  The  boundary  conditions 
used  to  achieve  this  are  described  in  Section  U,  together  with  the 
inital  conditions. 

Section  5 provides  a brief  description  of  the  numerical 
methods  which  we  have  employed  in  NRFL02  in  order  to  obtain  accurate  and 
convergent  results.  The  code  has  been  designed  to  be  economical  of 
machine  time,  despite  the  complexity  of  the  problem.  Implicit  time- 
stepping is  used  to  allow  larger  time  steps  as  we  evolve  the  solutions 
to  a steady  state.  A non-uniform  computational  mesh  is  used  to  achieve 
good  spatial  resolution  of  the  fairly  complicated  flow  structures. 
Second-order  accurate  finite-differences  are  used  throughout. 
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In  Section  6,  we  propose  and  discuss  a fairly  simple  labora- 
tory experimental  simulation  of  stratified  turbulence  driven  by  inflows 
and  outflows,  in  support  of  our  theoretical  studies.  The  experiment  is 
designed  to  investigate  recirculation  between  the  inlet  and  outlet 
ports,  a matter  of  critical  concern  in  this  effort. 

In  the  final  Section  7,  we  present  numerical  results  obtained 
by  applying  our  code  NRFL02  to  the  proposed  laboratory  experimental 
simulation, for  a range  of  flow  rates.  The  results  show  that  for  low 
flow  speeds,  the  inflow  and  outflow  drive  relatively  independent  fluid 
motions  which  are  practically  horizontal.  Doubling  the  flow  speed 
produces  significant  turbulent  recirculation  between  the  inflow  and 
outflow  ports,  and  doubling  it  again  produces  major  recirculation. 

Such  simplified  experiments  should  play  an  important  role  in  the  design 
of  actual  OTPP  prototypes. 
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2.  EQUATIONS  AND  TURBULENCE  MODEL 

In  this  section,  we  present  the  equations  used  in  our  NRFL02 
code  to  describe  these  stratified  turbulent  flows.  The  resulting 
equations  (13)  are  listed  at  the  end  of  the  section,  so  that  the 
discussion  can  be  omitted  on  a first  reading. 

Using  the  usual  Boussinesq  approximation  and  a standard 
notation,  the  full  equations  of  heat,  momentum,  and  continuity  are 
written  as : 

v • (i*> 

3tui  ■ + »j,i>  - Vl}.J  ' <lb) 

- o . (lc) 

Here  u^  is  the  velocity  vector,  T is  temperature,  and  P is  a modified 

pressure.  Further,  V and  K are  respectively  the  kinematic  viscosity 

and  the  thermal  diffusivity,  and  b is  the  upward  buoyancy  force 

gSp/p,  taken  as  gaT,  where  a,  the  coefficient  at  thermal  expansion,  is 

assumed  constant.  The  third  coordinate  (x^)  is  upwards,  the  usual 

summation  convention  is  implied,  and  all  spatial  derivatives  are 

indicated  by  the  suffix  notation,  so  that  u . denotes  3u  /3x. . 

J » i J i 

We  will  express  the  velocity  field  as: 


where  u^  is  the  average  of  an  ensemble  of  flow  realizations,  and  uT 
is  a particular  turbulent  fluctuation  from  this  average.  The  other 
variables  like  temperature  will  be  written  in  the  same  form.  Then 
taking  the  ensemble  average  of  equations  (l)  yields: 
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v -Hj-  v-^'Kj  • 

3tui  ■ bS13  + {-?Si3  + + “J>t>  - Vi 

'^Khi  • 


(2a) 


(2b) 


i.i 


(2c) 


Here,  the  long  overbars  denote  ensemble  averages  of  products  of  the 
turbulent  fluctuations  uT  and  T".  Thus,  -uJt7  is  the  turbulent  heat 
flux,  and  -uju'  is  the  Reynolds  stress.  Subtracting  equations  (2) 
from  equations  (l)  yields  the  turbulent  fluctuation  equations. 


a t“ 

t 


v;  ■ b'613  + {-p'su + "K.j + »;,!> 


Vi 


(3a) 


(3b) 


i.i 


(3c) 


It  is  not  practicable  to  solve  equations  (l)  or  (3)  for  a 
very  large  number  of  flow  realizations,  in  order  to  obtain  the  mean 
temperature  distribution  T and  the  mean  flow  ^ by  averaging.  Instead, 
statistical  turbulence  models  (Piacsek,  Toomre,  and  Roberts,  1975)  are 
used  to  obtain  equations  for  u^T7  and  which  are  partly  empirical, 
but  are  motivated  by  equations  (3). 
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In  sero-order  turbulence  modelling,  the  turbulent  heat  flux 
and  the  Reynolds  stress  are  determined  by  the . equations 

-u^T”1  = (K  - <)T  . , (U) 

J » J 

-vT  * (K  - * "i.i’  -7%  • (ltb) 

where  the  turbulent  diffusivity  K is  an  imposed  function  of  position, 
and  E is  the  turbulent  kinetic  energy  density  %uTu'.  In  general,  K can 
be  a tensor  rather  than  a scalar,  and  can  be  different  in  equations 
( Ua ) and  ( Ub ) . 

In  second-order  statistical  turbulence  modelling,  the 
quantities  u'T'  and  u'uf  obey  partial  differential  equations,  with  time 

1 J 1 

derivative  terms,  advection  terms,  generation  terms,  diffusion  terms, 
and  decay  terms.  These  equations  are  as  far  as  possible  derived  from 
equations  (3),  but  are  partly  empirical,  involving  concepts  such  as 
the  length  scale  or  turnover  time  of  the  turbulence,  as  in  the  following 
discussion  of  first-order  closure. 

We  adopted  for  the  code  NRFL02  a first-order  statistical 
closure  model  using  equations  (U),  with  K determined  from  an  imposed 
length  scale  L and  a turbulent  kinetic  energy  equation.  This  equation 
is  obtained  by  multiplying  equation  (3b)  by  ur,  taking  the  ensemble 
average,  and  modelling  certain  terms.  With 

E = \ uiui  ’ (5) 


the  resulting  turbulent  kinetic  energy  equation  is 
3tE  « ^ - ultj  + {-uj'(P>  E') 


(6) 


+ v ui(ui,j  +Uj.i)  ■ ujE },j  • iiuli  +ui,i)! 
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The  first  two  terms  on  the  right  of  equation  (6)  represent, 
respectively,  the  loss  of  turbulent  kinetic  energy  in  creating  mean 
potential  energy,  and  the  generation  of  turbulent  kinetic  energy  from 
the  shear  of  the  mean  flow.  We  write  the  buoyancy-related  term  u'b^ 
as  gotu^T'.  Equations  (4)  can  then  be  used  to  approximate  these  first 
two  terms.  The  third  term  can  be  interpreted  as  the  divergence  of  a 
flux  of  turbulent  kinetic  energy;  the  first  two  parts  are  modelled  as 

{ + + “jV  I .3  ■ Hj  }.s  • (7) 


The  last  term  in  equation  (6)  represents  the  loss  of  turbulent  kinetic 
energy  by  a cascade  through  the  spectrum  to  small  length  scales  where 
it  is  dissipated  by  viscosity.  This  process  occurs  on  a time  scale 
proportional  to  the  turnover  time  for  the  largest  eddies,  so  that  the 
last  term  is  modelled  as 


+ 


(3) 


where  c ^ is  a numerical  coefficient  and  t is  a turnover  time  scale. 

In  unstratified  turbulence  with  first-order  modelling, 

- •<  — L 

usually  K = LE  and  x = LE  , where  the  length  scale  L is  an  imposed 
constant  or  function  of  position,  or  is  determined  from  another  empi- 
rical equation.  When  dealing  with  turbulence  in  a stratified  medium, 
the  importance  of  the  stratification  is  measured  by  the  dimensionless 
ratio  N2L2/E  where  N2  is  b . and  N is  the  Brunt-Vaisala  frequency  of 

9 3 

internal  waves.  When  this  ratio  is  very  small,  the  turbulent  eddies 
turn  over  in  a time  short  compared  with  1/N,  and  the  stratification  is 
negligible.  When  the  ratio  is  very  large,  the  random  motions  are 
more  like  internal  waves  than  turbulent  eddies  and  the  system  of 
equations  (3)  is  effectively  linear  in  the  turbulent  fluctuations. 
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This  has  two  consequences . There  is  practically  no  nonlinear  cascade 
of  turbulent  kinetic  energy  to  small  length  scales  where  it  is  lost 
through  viscous  dissipation.  Also,  there  is  effectively  no  vertical 
diffusion  of  heat;  a vertically  displaced  fluid  element  bounces  back 
to  its  undisturbed  elevation  without  mixing.  When  the  dimensionless 
ratio  N2L/E  is  of  order  unity,  the  random  internal  waves  are  sufficiently 
nonlinear  to  break  at  irregular  intervals,  resulting  in  some  vertical 
diffusion  and  in  some  turbulent  kinetic  energy  loss  through  the  cascade. 

—h  — -h 

The  unstratified  formulae  K = LE^  and  T = LE  2 therefore 
give  too  large  a diffusivity  and  too  small  a time  scale  for  energy 
cascade,  when  stratification  is  important  end  the  dimensionless  ratio 
is  of  order  unity  or  larger.  In  addition,  when  this  K expression  is 
used  in  modelling,  the  first  term  on  the  right-hand,  the  turbulent 
kinetic  energy  density  E,  decays  locally  to  zero  in  a finite  time. 

For  this  reason,  we  have  adopted  for  the  code  NRFL02  the  following 
expressions  for  K and  T, 

K = LE^Cl  + c N2L2/E)_1  , (9) 

s 

x = L2/K  , (10) 

with  L an  imposed  constant  length  scale  and  cg  a constant. 

For  our  modelling  of  the  external  flow  for  an  OTPP,  we  have 
made  further  modifications  to  the  system  of  equations  (2)  and  (i+)  to 
(10)  to  give  steady  solutions  when  the  dimensionless  ratio  N2L2/E  is 
large.  The  mean  flow  variables  T,  u^,  and  E must  be  steady,  from  their 
definition  as  the  average  over  an  ensemble  of  flows  maintained  by 
steady  boundary  conditions.  And  in  cases  with  reasonably  strong  mean  flow 
and  turbulence,  we  have  obtained  steady  solutions  by  time-stepping  the 
above  equations.  But  in  other  cases,  we  obtained  a finite  amplitude 
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oscillation;  it  appears  that  the  unmodified  syst  m of  equations 
does  have  a steady  solution,  but  this  solution  is  unstable  to  small 
disturbances  when  the  turbulent  diffusivity  K is  small. 

We  have  considered  three  ways  of  resolving  this  difficulty 
that  the  mean  flow  variables  oscillate  in  time  when  the  dimensionless 
ratio  N2L2/E  is  small.  The  first  is  the  rational  approach  of  time 
averaging  the  oscillations , but  such  a procedure  is  of  questionable 
accuracy.  The  second  approach  is  to  look  for  steady  solutions  to  the 
same  equations  by  a method  which  does  not  involve  marching  them 
forward  in  time,  thus  avoiding  the  instability.  The  third  method  is 
to  modify  the  equations  and  the  turbulence  modelling  in  order  to 
suppress  the  instability.  This  could  be  achieved  by  a large  and  totally 
artificial  increase  in  the  turbulent  diffusivity  K;  the  results  will 
then  be  steady,  but  numerically  incorrect.  A more  subtle  change  in 
the  equations  is  required,  and  it  is  this  approach  which  we  have 
implemented  in  our  NRFL02  code. 

To  achieve  this , we  have  added  the  term 


-(c  KN2/E)u, 


(11) 


on  the  right-hand  side  of  the  u^  equation.  The  term  taken  alone 

would  make  u,  tend  to  zero  on  a time  scale  of  (c  KN2/E)  1 . When 
J w 

used  in  combination  with  the  other  terms , and  with  the  consistent 
addition  of  the  term 


(c  KN2/E)u.2 
w 3 


(12) 


to  the  right-hand  side  of  the  turbulent  kinetic  energy  equation  (6),  it 
serves  to  suppress  the  fluctuating  vertical  motions  associated  with  the 
instability  of  the  mean  flow,  and  to  generate  turbulent  kinetic  energy 
instead.  With  these  added  terms  (11)  and  (12),  we  have  been  able  to 
obtain  steady  mean  flows  for  a wide  range  of  OTPP  flow  simulations. 
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The  resulting  equations  used  in  NRFL02  are 


V ' Ki  - V >.j  ' (Ua) 

3tui  = (gaT  " cwo2u3/E)6i3  + {"Vij 

+ C<“i.J  +“3.1)  ‘ Vil.i  1 (l3b) 

u =0  , (13c) 


9tE 


KfNi  + j(ui,j + uj,i)!-cfE/L! 

+ VS*'*}*  Rj  - y}  , 


(13d) 


N2  = gaT  3 , (13e) 

K = LE^/d  + csN2L2/E)  , (13f ) 

where  P_  is  P + |e  in  equation  (13b),  and  a term  ~|v(u  + u ,)2  has 

3 c 1 » J J > 1 

been  dropped  from  equation  (13d)  on  the  assumption  v«K. 

Equations  (13)  involve  four  parameters,  the  length  scale  L 

and  the  dimensionless  constants  cc,  c and  c . These  must  be  deter- 

f s w 

mined  by  tuning  their  values  so  that  the  numerical  solutions  are  in 
satisfactory  agreement  with  the  results  of  laboratory  experiments. 

The  thermal  expansion  coefficient  a is  taken  as  1.3  x 10  4/°F,  the 
appropriate  value  for  temperatures  near  ^0°F. 
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3. 


HE  OCEAN  THERMAL  POWER  PLANT  MODEL 


In  the  code  NRFL02,  the  ocean  thermal  power  plant  and  its 
inflows  and  outflows  are  modelled  by  boundary  conditions  on  the  left 
boundary  of  the  rectangular  domain  - D<z<0,  0<x<YD.  Note  that  in 
describing  our  two-dimensional  application  of  equations  (13)  we  use  x 
as  the  horizontal  coordinate  and  z as  the  vertical  coordinate,  in 
place  of  x^  and  x^,  and  u and  w for  u^  and  u^.  None  of  the  currently 
suggested  plant  designs  (Dugger,  1975)  are  described  exactly  by  such 
a two-dimensional  model,  but  it  does  serve  to  evaluate  their  primary 
flow  characteristics. 


For  each  of  four  inflows  and  outflows,  the  code  uses  a 
center  height  z^,  a radius  r^>  a signed  flux  amplitude  a^,  an  angle 
above  the  horizontal  plane,  a turbulence  constant  E^,  and  a temper- 
ature T^,  together  with  the  function 


f(x)  = 


for  I x | < 1 
for  j x | > 1 


(14) 


The  value  of  u on  the  left  boundary  at  height  z is 


4 


Thus  u is  zero  except  close  to  one  of  the  z^  values,  and  is  negative 
when  a^  is  negative  (a  plant  inflow)  and  positive  when  a^  is  positive 

(a  plant  outflow).  The  function  f (x)  has  been  chosen  to  make  the 

— 2 2 
boundary  value  of  u a smooth  function  of  z;  replacement  of  (1-x  ) 

by  1 in  equation  (14)  may  give  a more  faithful  description  of  certain 

plant  designs  or  experimental  simulations.  If  fewer  than  four  plant 

inflows  and  outflows  are  to  be  modelled,  then  some  of  the  a^  values 

can  be  set  to  zero. 
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For  a plant  outflow,  where  fluid  is  entering  the  computational 
region,  boundary  values  must  be  imposed  for  w,  T,  and  E.  The  values 
are 


w(z) 

* u tan  a. 

k 

9 

(16a) 

E(z) 

- V*V2 

9 

(16b) 

T(z) 

■ \ • 

(16c) 

Occasionally,  slightly  more  complicated  formulae  are  used  in  NRFL02  to 
allow  for  the  possibility  of  overlapping  outflow  regions. 

Since  the  power  plant  uses  warm  surface  water  to  vaporise  a 
fixed  flux  of  working  fluid,  the  warm  outflow  temperature  is  lower  than 
the  mean  warm  inflow  temperature  by  a fixed  amount.  Similarly,  the  cold 
water  inflow  is  used  to  condense  the  working  fluid,  and  its  temperature 
is  increased  a fixed  amount  before  it  is  pumped  out  as  the  cold  out- 
flow. The  code  therefore  evaluates  the  mean  temperature  T^  for  each  plant 
inflow  region  (a^  negative), 


and  adds  or  subtracts  a given  temperature  increment  to  determine  the 
corresponding  outflow  temperature  T^.  This  is  optional;  the  code  can 
also  use  fixed  T^  values  for  plant  outflows  as  in  the  experimental 
simulation  described  in  Section  6. 

We  point  out  here  that  the  NRFL02  code  has  the  flexibility  to 
describe  a variety  of  plant  designs.  Further,  it  can  be  used  to  eval- 
uate a portion  of  a given  design,  to  evaluate  for  example  just  the  cold 
inflow,  or  just  the  warm  inflow  and  the  two  outflows.  This  flexibility 
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in  the  code  allows  it  to  be  used  to  describe  other  configurations,  such 
as  the  heated  outfall  from  a conventional  power  plant,  or  the  simple 
experimental  simulation  presented  in  Sections  6 and  7. 


h. 


BOUNDARY  AND  INITIAL  CONDITIONS 


In  many  hydrodynamic  simulations,  the  boundary  conditions  of 
symmetry  or  periodicity  are  essentially  passive,  and  are  only  minor 
aspects  of  the  problem.  In  our  near-field  simulations  for  ocean 
thermal  power  plants,  it  is  the  boundary  conditions  on  the  left 
boundary,  described  in  Section  3,  which  drive  the  flow,  and  the  remain- 
ing boundary  conditions  which  control  it. 

At  large  distances  from  an  OTPP  operating  in  a stratified 
ocean,  th«!  flow  is  essentially  horizontal,  and  the  outflow  water  from 
the  near-field  region  has  found  its  own  density  level.  Thus  a primary 
requirement  for  a numerical  simulation  of  the  near  flow  is  a set  of 

boundary  conditions  at  the  right-hand  edge  of  a finite  computational 
domain,  which  make  this  boundary  permeable  and  allow  horizontal  inflow 
and  outflow  a*  the  correct  density  levels. 

A second  requirement  is  that  the  boundary  conditions  where 
the  flow  is  out  of  the  computational  region  should  be  passive.  They 
should  not  produce  significant  outflow  boundary  layers,  nor  should 
failure  to  resolve  these  boundary  layers  numerically  lead  to  upstream 
influence. 

A third  desirable  feature  for  these  boundary  conditions  is 
that  during  time-stepping  towards  a steady  solution  of  the  equations, 
they  should  make  the  boundary  permeable  to  internal  waves,  and  not 
allow  them  to  reflect.  As  stated  in  Section  2,  there  are  difficulties 
in  these  computations  over  instability  of  the  mean  flow  solution;  a 
boundary  condition  which  reflects  oscillations  will  increase  these 
difficulties. 

In  our  code  NRFL02,  we  assume  the  ambient  ocean  T,  w,  and  E 
distributions  are 

Ta(z)  - Tt  + Tr  tan-‘j(z  + d^/zj  , (18a) 

wfl  - 0 , (18b) 
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E (z)  - E exp  (z/z  ) 

3 0 C 


(18c) 


i 

f 


where  Tc  is  the  thermocline  temperature,  T the  thermocline  temperature 
range,  z£  the  thermocline  thickness,  and  d£  the  thermocline  depth,  while 
the  turbulent  kinetic  energy  density  amplitude  and  depth  scale  are  Eq 
and  Zg.  respectively. 

Equation  (18a)  implies  a temperature  variation  from 
Tt  - ttT^/2  at  great  depths  to  Tt  + ttT^/2  at  large  z (though  z is  of 
course  confined  to  negative  values).  Parameter  values  can  be  chosen 
so  that  this  profile  is  a satisfactory  approximation  to  most  real 
temperature  profiles.  According  to  equation  (13e) , the  corresponding 
Vaisala  frequency  N of  internal  gravity  waves  in  the  undisturbed  ocean 
is  given  by 

N2  - gaTrztjz2  + (z  + dt)2j  . (19) 


Equation  (18c)  is  a reasonable  choice,  since  the  ambient  turbulence 
in  the  upper  layers  of  the  ocean  is  mostly  due  to  the  wind-driven 
surface  waves,  and  decreases  with  depth.  In  fact,  the  rate  of  decrease 
is  more  rapid  when  N2  is  large.  But  since  the  ambient  turbulence  is 
much  smaller  than  the  turbulence  produced  by  the  plant,  our  solutions 
are  relatively  insensitive  to  Eq  and  z^.,  and  thus  equation  (18c)  is 
adequate. 

At  the  ocean  surface  z ■*  0 and  at  the  lower  boundary  of  the 
computational  domain  z - -D,  we  apply  the  boundary  conditions 


w ■ 0 


(19a) 


E 


z 


= 0 


» 


where  u 

.z 


denotes  3u/3z. 
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(19b) 

(19c) 

(I9d) 


Thus  there  is  zero  normal  flow,  and  zero  turbulent  flux  of  horizontal 
momentum,  heat,  and  turbulent  kinetic  energy  across  these  boundaries. 

In  applying  equations  (19b),  (19c)  and  (19d)  at  the  surface,  we  are 
neglecting  respectively  the  momentum  flux  across  the  surface  produced 
by  wind  stress,  the  heat  flux  (due  to  solar  heating,  sensible  heat 
transfer,  and  evaporative  and  radiative  cooling),  and  the  generation 
of  a flux  of  turbulence  at  the  surface  by  wind  and  waves.  Each 
approximation  is  justified  by  the  fact  that  the  order  of  magnitude  of 
each  of  these  fluxes  in  our  near-field  computation  is  much  larger  than 
the  geophysical  surface  fluxes. 

We  next  describe  our  boundary  conditions  at  the  left  and  right 
boundaries  of  the  computational  domain.  On  the  left  boundary,  the 
normal  velocity  u is  given  by  equation  (15) , representing  the  OTPP 
inflows  and  outflows.  In  order  to  satisfy  our  first  requirement 
that  the  right  boundary  be  permeable,  allowing  horizontal  inflow  and 
outflow  with  the  appropriate  temperature  for  each  depth,  we  have 
applied  the  equation 

- 8CTa  , (20) 

at  x * yD,  where  T is  the  ambient  temperature  distribution  (18a). 

This  implies  a pressure  distribution  in  hydrostatic  balance  with  the 
buoyancy  force  associated  with  the  ambient  temperature  distribution. 
From  equation  (13b),  this  balance  must  apply  at  large  distances  from 
the  plant,  when  the  flow  is  horizontal  and  the  turbulent  kinetic  energy 
density  and  the  turbulent  diffuslvity  are  negligibly  small.  With  the 
boundary  condition  (20),  the  only  possible  horizontal  flow  solution 
near  the  boundary  x * yD  must  have  the  mean  temperature  distribution 
T(x,z)  - Ta(z). 

We  further  need  boundary  conditions  on  T,  w,  and  E at  the 
side  boundaries.  These  boundary  conditions  take  a different  form 
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depending  on  whether  the  horizontal  flow  is  into  or  out  of  the  compu- 
tational domain.  For  outward  flow,  the  boundary  conditions  must  be 
passive,  to  satisfy  our  requirement  that  no  significant  outflow  boundary 
layer  be  formed.  For  both  outward  and  zero  normal  flow,  we  impose 
zero  normal  derivative  on  both  vertical  boundaries.  Thus  for  u not 
inward. 


T - 0 

,x 


,x 


»X 


0 


(21a) 

(21b) 

(21c) 


on  x * 0 and  x » yD.  For  inward  flow,  the  values  of  T,  w,  and  E are 
determined  by  equations  (16)  for  the  left  boundary  (the  OTPP  outflow 
values),  and  by  equations  (18)  for  the  right  boundary  (the  ambient 
ocean  values) . 

Our  third  requirement,  that  the  right  boundary  be  permeable 
to  internal  waves,  is  not  met  by  equation  (20)  and  the  other  boundary 
conditions  on  the  right.  We  have  therefore  added  the  further  term 


-a(x, z)w 


(22a) 


to  the  w equation.  Here  a is  only  non-zero  close  to  the  right  hand 
boundary,  where  it  serves  to  damp  vertical  motions.  Incoming  internal 
waves  are  absorbed  by  this  "porosity"  distribution  (Piacsek  and 
Roberts,  1975).  For  o,  we  have  used  the  equation 


Cf(x,z) 


jN  gj  8(x/yD 


1) 


(22b) 


where 


g(0  " exp  J -sinh(sinh  C)  j 


(22c) 
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The  initial  conditions  for  equations  ( 1 3 ) can  be  stated 
readily.  The  NRFL02  calculations  are  initiated  by  asserting  ambient 
ocean  conditions  (l8)  throughout  the  computational  domain.  Further, 
the  initial  horizontal  velocity  profile  at  all  values  of  x is  the  same 
function  of  z as  that  prescribed  on  the  left  boundary  by  equation  (15). 
We  have  chosen  these  initial  conditions,  with  horizontal  flow  and  very 
weak  turbulence,  so  that  any  recirculation  which  occurs  is  not  due  to 
the  initial  conditions.  In  fact,  our  calculations  have  not  provided 
evidence  for  any  hysteresis  effects  of  the  initial  conditions  on  the 
final  steady  solution. 
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5. 


NUMERICAL  PROCEDURE  FOR  SOLVING  THE  EQUATIONS 


In  this  section  we  give  a brief  summary  of  the  numerical 
methods  we  have  used  to  solve  equations  (13) , with  the  initial  and 
boundary  conditions  described  in  Sections  3 and  4.  Since  prolonged 
calculations  are  envisaged  for  a large  number  of  flow  situations 
and  for  different  parameter  ranges,  the  numerical  procedures  must  be 
very  efficient.  Further,  it  is  desirable  that  the  computer  code  used 
to  implement  these  methods  should  possess  a structure  with  a high 
degree  of  clarity  and  flexibility. 


We  begin  by  describing  the  non-uniform  computational  mesh 
used  in  our  finite-difference  representation  of  the  equations.  We  use 
a small  mesh  spacing  only  near  the  left  boundary  and  near  the  plant 
inflows  and  outflows,  where  it  is  needed  to  resolve  the  small  length 
scales  present.  Elsewhere  we  use  a larger  mesh  interval,  to  economize 
on  mesh  points.  The  mesh  spacing  varies  smoothly. 


Our  computational  mesh  is  staggered,  as  shown  in  Figure  1, 

for  greater  efficiency  in  representing  equations  (13).  The  variables 

T,  E,  P-  and  K are  defined  at  the  mesh  points  (x  , z ) for  integer 

m n 

values  of  m and  n,  while  u is  defined  for  integer  n and  half-odd  m, 

and  w for  integer  m and  half-odd  n.  It  has  long  been  known  that  this 

mesh  is  much  more  convenient  and  accurate  for  calculations  involving 

the  primitive  variables.  Most  of  the  required  first  derivatives  can 

be  represented  using  neighboring  mesh  points,  thus  u + w in 

» x ,2  

equation  (13c)  is  conveniently  represented  at  the  T points,  and  P 

___  L yX 

and  Pr  are  conveniently  represented  at  the  u and  w points. 

£• , Z 


Our  boundaries  are  placed  between  the  mesh  lines,  rather  than 
on  them,  in  representing  the  boundary  conditions.  As  shown  in  Figure 
1,  the  side  boundaries  are  at  and  x^  and  the  bottom  and  surface 
boundaries  are  at  z ^ and  z^_^.  This  choice  puts  the  normal  velocity 
components  on  the  boundary,  and  makes  for  high  accuracy  in  representing 
the  normal  derivatives  of  the  other  variables. 


Next  we  describe  the  numerical  representation  used  for  the 
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Figure  1.  The  non-uniform  staggered  computational  mesh.  T,  E,  P and  K 
are  represented  by  their  values  at  the  T points,  while  u and 
w are  represented  by  their  values  at  the  staggered  mesh 
points  indicated.  The  boundaries  are  at  x^,  x^jj,  zju, 
zn-*5'  This  leads  to  greater  accuracy,  convenience,  ana 
efficiency  in  representing  the  equations  and  boundary  condi- 


tions. 


spatial  derivatives  in  equations  (13).  For  every  term  but  one,  we 
employ  the  simplest  second-order  accurate  central  difference  repre- 
sentation available,  with  averaging  from  the  neighboring  two  or  four 
mesh  points  as  required.  The  resulting  finite  difference  equations 
conserve  the  volume  flux,  heat,  and  momentum,  as  implied  by  equations 
(13).  The  one  deviation  from  central  differences  is  the  use  of  upwind 
differencing  in  representing  the  advection  term  in  the  turbulent 
kinetic  energy  equations  (13d).  For  this  term 


(u . E)  . = uE  + wE 
1 . j ,x  ,z 


we  use  the  representation 


(uE  ) = j(u  + | u |)(E  - E ,)/6x  i 

, x m * 1 m m- 1 m-* 


+ -(u  - j u j ) (E_ 


E )/&c 


where  the  suffix  n for  the  z mesh  is  understood  and  u denotes 
jju^  ^ + u^^  ( . This  representation  is  only  first-order  accurate; 
it  has  been  adopted  to  avoid  negative  E values  where  there  is  flow 
from  a region  of  small  £ to  a region  of  large  E and  the  intervening 
layer,  which  can  be  very  thin,  is  not  sufficiently  resolved. 

The  boundary  conditions  described  in  Sections  3 and  4 are 
similarly  represented  using  the  simplest  central  difference  method, 
with  one  exception.  Where  the  flow  is  out  of  the  computational 
domain  (to  the  left  or  the  right),  the  T,  w and  E equations  allow  an 
outflow  boundary  layer  of  thickness  K/|u|.  When  this  is  smaller  than 
6x/4,  our  resolution  is  inadequate,  and  there  is  a tendency  for  mesh 
separation,  with  the  variables  at  even  mesh  points  separating  from  the 
variables  at  odd  mesh  points.  This  tendency  is  eliminated,  and  results 
of  greater  accuracy  are  obtained,  by  replacing  the  boundary  condition 
of  zero  first  derivative  with  zero  second  normal  derivative,  even 
when  the  representation  is  only  first  order. 
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We  have  chosen  to  represent  the  time  evolution  of  equations 
(13)  with  second  order  accuracy,  in  spite  of  the  fact  that  we  are 
mainly  concerned  with  time-stepping  to  a steady  solution.  Due  to  the 
extreme  nonlinearity  of  the  equations  and  our  desire  to  avoid  the 
artificial  insta'oility  or  dissipation  associated  with  other  explicit 
representations,  we  were  forced  to  use  a "leapfrog"  method  for  time- 
stepping equations  (13).  The  diffusion  terms  were  treated  by  an 
alternating  direction  implicit  method.  We  write  the  equations  (l3a,b,d) 
as : 

3 f = Q + Xf  + Zf  , (25) 

where  f can  be  T,  u,  w,  cr  E; 

X and  Z are  the  diffusion  terms  involving  two  differentiations 
with  respect  to  x and  z; 

Q represents  the  rest  of  the  terms  in  the  equation. 

Then  for  time  kSt , where  k is  an  integer,  our  representation  of 
equation  (25)  is 


(f*  - fk_1)/6t  = Qk  + Xkfk_1  + Zkf* 

» 

(26a) 

(fk+1  - f*)/6t  = Qk  + Xkfk+1  + Zkf* 

(26b) 

At  the  beginning  of  the  time  step  all  quantities  f and  f are 
known,  and  thus  Q can  be  evaluated.  The  z-dif fusion  is  done  implicitly 

in  the  first  stage,  and  the  x-dif fusion  in  the  second.  The  intermediate 

* k 

quantity  f is  a second-order  approximation  to  f . 

Equations  (26)  are  subject  to  the  well-known  leapfrog 

instability  when  applied  to  the  turbulent  kinetic  energy  density 

equation  (13d),  because  the  N2  and  c,E/L2  terms  in  the  first  bracket 

are  decay  terms.  To  avoid  this  instability,  in  which  E*  for  even  k 

separates  from  Ek  for  odd  k,  we  subtract  the  term  2ak(E*  - ’Ek)  from 

the  right  hand  side  of  equation  (26a) , and  the  corresponding  term 
k — ~k"^  ^ j 

a (E  - 2E  + E ) from  the  right  hand  side  of  equation  (26b) . 
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Here 


2ak+  ]-N2  + i(u,  , + u.  )*  - c -E/L2  + c N2w2/  E i 
( 1 > J J > i I w / 


3K/3E 


-K(c ,/L2  + c N2w2/E2)  >0 
I w 


(27) 


A suitable  a expression  is  evaluated  and  applied  in  the  code. 

— k 

The  pressure  P£  remains  to  be  determined,  in  order  to  apply 
equations  (26)  to  equation  (13b).  We  write 


(28) 


where  is  obtained  by  extrapolation  from  previous  values, 

■ 2*T  - pe~2  • 


(29) 


so  that  Pk  remains  unknown,  but  is  of  order  6 1 2 . Then  Pk  is  used  in 

evaluating  Q in  equation  (25),  and  equations  (26)  determine  u and  w- 

— k+i  — k+J  ^ ^ 

which  are  third  order  approximations  to  u and  w , and  do  not  satisfy 

1c  k 

equation  (13c).  We  then  add  the  effects  of  P„  _ and  P , applied 
for  a time  interval  26t, 


(uk+1  - V/26t  " - P2,x 


(30a) 


(wk+i-  w9)/26t  = 


- p: 


Z,z 


(30b) 


Taking  the  divergence  and  using  equation  (13c)  gives 


,IP2  ' <“2,x  + »2,x)/26t 


(31) 


This  Poisson  equation  is  solved  by  a generalized  form  of  Wachspress 
alternating-direction  implicit  iteration  (Varga,  1962,  Roberts  and 
Piacsek,  1975)  using  four  iterations,  while  imposing  the  boundary 
conditions , 
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lr 

= 0,  on  x = YD, 


(32a) 


1, 

SP^/Sn  = 0,  on  the  other  boundaries.  (32b) 

The  resulting  P solution  is  substituted  in  equations  (28)  and  (30)  to 

~k  r 1 -i 

determine  P£  (_for  application  of  equation  (29)J  and  the  new  velocity 

components . 

With  the  numerical  method  described  by  equations  (26),  our 
code  is  still  subject  to  the  stability  requirement 


(N  + |u|/<5x  + |w|/6z)6t  <1 


(33) 


Since  6t  is  the  same  at  every  mesh  point  and  time  step,  an  estimate 
must  be  made  at  the  beginning  of  the  calculation  of  the  maximum  value 
of  the  bracket  which  is  likely  to  be  encountered.  The  code  uses 


6t  = c(N 


max 


u /6x  . ) 
max  min 


-1 


(34) 


where  N and  u are  the  maxima  of  the  ambient  stratification  (19) 
max  max 

and  the  left  boundary  velocity  (15)  respectively,  and  c is  an  input 
constant . 

Leanfrog  time-stepping  (26)  is  subject  to  a weak  computational 
instability,  with  the  values  separating  at  even  and  odd  time-steps,  even 
in  the  absence  of  decay  terms.  These  parasite  solutions  are  avoided  by 
re-iritializing  the  calculation  every  few  time-steps.  For  initialization, 

in  terms  of  equation  (26),  any  previous  f is  forgotten,  and  f is  used 

k rk  k k+1 

to  evaluate  Q , a , and  Z in  order  to  obtain  an  approximate  f^  , with 

second-order  accuracy.  Then  a new  calculation  can  be  started  off  from 
„k-l 

f using 
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6. 


A PROPOSED  LABORATORY  SIMULATION 


Our  code  NRFL02  calculates  the  flow  in  a rectangular  domain 
representing  the  near-flow  region  closest  to  an  ocean  thermal  power 
plant,  where  the  flow  deviates  from  the  horizontal  and  where  the 
turbulence  generated  by  the  plant  is  important.  The  rectangular  region 
has  impermeable  boundaries  at  the  top  and  bottom.  The  horizontal 
velocity  is  given  on  the  left,  with  the  temperature,  vertical  velocity, 
and  turbulent  energy  density  also  being  imposed  where  this  velocity  is 
out  of  the  plant  and  into  the  region.  The  passive  boundary  conditions 
on  the  right  allow  horizontal  inflow  from  and  outflow  to  a large  fluid 
region  with  assumed  temperature  and  turbulent  kinetic  energy  profiles. 

This  code  NRFL02  has  four  parameters : the  turbulence  length 
scale  L,  and  the  three  dimensionless  parameters  c^,  cg , and  c^.  Before 
the  code  can  be  used  with  any  confidence  to  simulate  the  external  flow 
of  an  ocean  thermal  power  plant , we  must  tune  the  values  of  these 
parameters  so  that  the  numerical  results  are  in  satisfactory  agreement 
with  the  results  of  laboratory  experiments.  Such  experiments  should  be, 
two-dimensional,  and  involve  the  simultaneous  effects  of  density 
stratification  and  of  turbulent  jet  inflows  and  outflows  from  the  region. 
To  our  knowledge , suitable  experiments  have  not  yet  been  done , but  two 
experimental  programs  recently  funded  by  ERDA  should  begin  to  correct 
this  deficiency. 

We  now  suggest  an  experiment  which  meets  these  conditions  and 
which  would  provide  a laboratory  examination  of  the  possibility  of 
near- flow  recirculation  of  the  outflow  water  in  OTPP  operation,  as  well 
as  of  the  far-field  environmental  effect.  In  addition,  this  experiment 
would  provide  data  for  verifying  and  tuning  our  code  NRFL02  and  its 
successors. 

The  proposed  experimental  configuration  for  this  OTPP  near- 
field simulation  is  outlined  in  Figure  2.  A long  rectangular  tank  is 
filled  with  water  to  a depth  of  U ft,  with  the  temperature  increasing 
linearly  from  i+0°F  at  the  bottom  to  80°F  at  the  surface.  At  the 
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FIGURE  2 


The  proposed  laboratory  experiment  configuration. 

A long  tank  is  initially  filled  with  water  to  a 
depth  of  1*  ft,  with  the  temperature  increasing  linearly 
from  U0°F  at  the  bottom  to  80°F  at  the  free  surface. 

Water  at  **0°F  is  pumped  steadily  in  through  the  bottom 
gate  on  the  left  wall,  and  an  equal  volume  flux  is 
removed  at  the  top.  The  inflow  drives  the  indicated 
turbulent  Jet,  entraining  water  from  above.  Further 
along  the  tank  (beyond  about  8 ft)  the  turbulence  becomes 
very  weak,  and  the  flow  is  essentially  horizontal,  as 
indicated.  The  outflow  from  the  tank  comes  partly  from 
this  far-field  horizontal  motion,  and  partly  from  the 
cold  inflow;  the  mean  outflow  temperature  can,  therefore, 
range  from  almost  80°F  for  very  slow  flow  to  U0°F  for  very 
fast  flow.  The  inflow  and  outflow  maximum  speeds 
U = 2.68  ft/sec,  1.3**  ft/sec,  and  0.67  ft/sec  are  studied 
numerically  in  Section  7. 


29 


beginning  of  the  experiment  the  water  is  at  rest.  Cold  water  at  l*0°F 
is  then  pumped  in  steadily  at  the  end  of  the  tank,  through  a gate  1 ft 
high  at  the  bottom,  and  taking  up  the  full  width  of  the  tank.  This 
flow  simulates  an  OTPP  outflow.  To  simulate  a plant  inflow,  the  same 
volume  flux  of  water  is  simultaneously  removed  from  the  same  end  of 
the  tank,  through  a similar  gate  extending  1 ft  down  from  the  surface. 
The  maximum  inflow  speed  U for  the  proposed  experiment  should  be  in 
the  range  from  0.2  ft/sec  to  5 ft/sec.  This  speed  range  brackets  the 
critical  value  described  in  the  simple  theoretical  discussion  of  the 
resulting  stratified  turbulent  flow  provided  in  the  remainder  of  this 
section. 

In  the  absence  of  density  stratification,  a turbulent  flow 
field  is  established  fairly  quickly  in  which  all  of  the  cold  inflow 
water  is  sucked  out  again  at  the  top  of  the  tank.  The  mean  flow  and 
the  turbulence  intensity  E become  very  small  at  large  distances  (more 
than  twice  the  depth)  from  the  end  of  the  tank.  Thus,  the  far-field 
flow  is  zero. 

With  a density  stratification,  the  situation  is  quite 
different.  At  large  distances  from  the  end  of  the  tank,  the  turbulence 
is  still  very  weak,  but  there  can  be  substantial  horizontal  motion. 

This  horizontal  flow  is  established  by  very  small  horizontal  density 
gradients,  and  is  limited  only  by  continuity  (13c)  and  by  the  stability 
requirement, 

N2  > (3u/9z)2  , (35) 

which  follows  from  the  first  two  terms  in  the  turbulent  kinetic  energy 
density  equation  (13d).  Physically,  when  this  equation  is  satisfied, 
more  turbulent  energy  is  lost  to  creating  potential  energy  than  is 
gained  from  the  shear.  Thus,  with  N constant  (from  equation  (13e),  N 
in  these  experiments  is  0.20U  sec  ^),  the  maximum  stable  far-field  u 
distribution  is 
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u 


N(z  + 2)  (x/X  - 1) 


(36) 


where  X is  the  length  of  the  tank,  and  the  total  flux  J udz  along 

-k 

the  tank  is  zero.  The  continuity  equation  (13c)  implies  a very  weak 
upwelling  velocity,  independent  of  x,  and  given  by  the  equation 

w = -Nz ( z + U)/2X  , (37) 

for  the  maximum  possible  u distribution  (36).  This  upwelling,  or  the 
distribution  corresponding  to  a u smaller  than  (36),  will  eventually 
modify  the  ambient  temperature  distribution  in  the  tank,  in  a time  of 
order  l/w. 

For  small  values  of  U,  therefore,  the  cold  inflow  Jet  becomes 
a horizontal  flow  as  the  turbulence  decays.  The  surface  tank  outflow 
simultaneously  draws  water  steadily  away  from  the  surface  layer.  The 
resulting  far-field  flow  is  smaller  than  that  given  by  equation  (36), 
and  is,  therefore,  stable.  The  average  outflow  temperature  is  almost 
80°F,  since  all  the  outflow  comes  from  near  the  surface  and  since 
minimal  turbulent  heat  transfer  occurs  between  the  warm  upper  layers 
and  the  cold  lower  layers . 

For  larger  U values,  in  a critical  range  and  above,  the  far- 
field  flow  will  be  approximately  given  by  equation  (36).  However,  this 
flow  is  too  small  to  carry  away  the  cold  inflow  and  to  provide  the  water 
for  the  surface  outflow.  The  remainder  of  the  cold  inflow  water  is 
sucked  up  and  out  with  the  surface  outflow.  (If  carried  over  to  an 
OTPP  design,  this  would  imply  recirculation  of  the  cold  plant  outflow 
into  the  warm  surface  intake.)  The  average  surface  outflow  temperature 
(17)  decreases  to  kO°F  as  U becomes  increasingly  large,  and  the  far- 
field  warm  surface  water  flow,  with  u negative  according  to  equation  (36), 
is  diluted  with  an  increasing  flux  of  cold  inflow  water. 
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According  to  equation  (36),  the  maximum  possible  far- field 

outward  cold  flux  J udz  is  2N,  for  x « X.  Now  the  cold  inflow  flux 

3 _ -U 

udz  is  XU,  where  U is  the  maximum  speed  and  X depends  on  the 

inflow  profile,  with  h < X < 1.  So  the  critical  U value,  for  which 
all  the  cold  inflow  can  flow  out  to  infinity,  is  2N/X.  This  is  about 
0.77  ft/sec  for  the  value  X = 8/15  used  in  Section  7- 

For  the  laboratory  experiment,  the  tank  width  should  be  at 
least  U ft,  to  reduce  three-dimensional  effects  on  the  mean  flow.  The 
length  X should  be  at  least  ho  ft,  to  allow  a reasonable  period  (about 
3x/2  sec)  for  establishment  of  a statistically  steady  turbulent  flow 
before  the  upwelling  (37)  significantly  modifies  the  temperature  and 
density  stratification.  The  temperature  and  the  velocity  components 
should  be  recorded,  as  functions  of  time,  at  a number  of  points  in  the 
flow,  in  order  to  evaluate  the  ensemble  average  (or  time  average) 
quantities  T,  u,  w,  u T , w T',  u'u',  u w , and  v'v'.  The  measurement 
points  should  be  mainly  in  the  first  8 ft,  but  should  include  a few 
points  further  along  the  tank  to  confirm  our  predictions  of  horizontal 
flow  and  small  turbulent  fluctuations.  Most  of  the  measurement  points 
should  be  in  the  x-z  plane  at  the  center  of  the  tank,  to  reduce  effects 
of  the  sidewalls.  These  sidewall  effects  can  be  estimated  by  adding 
measurement  points  across  the  whole  width  of  the  tank  at  selected 
values  of  x and  z. 
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7. 


NUMERICAL  RESULTS  FOR  THE  PROPOSED 
EXPERIMENTAL  SIMULATION 


We  have  applied  our  code  NRFL02  to  the  proposed  experimental 
simulation  of  stratified  turbulence  driven  by  fluid  jets  and  sinks,  as 
described  in  Section  6.  This  numerical  simulation  of  a proposed  experi- 
ment has  served  three  purposes.  First,  it  has  provided  a simple  and 
economical  problem  for  code  development,  and  for  improving  and  optimizing 
the  numerical  methods  of  Section  5,  since  adequate  resolution  and 
convergence  can  be  obtained  with  very  few  mesh  points  (about  12x12)  and 
time-steps  (about  100).  Secondly,  the  numerical  simulation  has  enabled 
us  to  study  the  effects  of  varying  the  dimensionless  parameters  cf,  c , 

and  c , and  the  ratio  of  the  turbulence  length  scale  L to  the  given 
w 

dimensions.  Such  studies  will  be  of  great  value  in  tuning  these 
parameters  to  give  results  agreeing  with  observations.  Thirdly,  the 
numerical  results  have  demonstrated  the  existence  of  the  far-field  flow 
maximum  (36)  and  the  corresponding  critical  speed  U for  inflow  to  the 
experimental  tank.  This  critical  speed  corresponds  to  an  upper  limit 
for  the  outflow  of  an  ocean  thermal  power  plant  prototype,  to  avoid 
recirculation  problems  and  thermal  efficiency  loss. 

For  our  NRFL02  simulation  of  the  proposed  experiment,  we 
took  the  horizontal  flow  u on  the  left  boundary  x=0  as 

u = u|f(z  + U)  - f(z)  J , (38) 

in  conformity  with  equations  (lM  and  (15).  Thus,  the  inlet  and  outlet 
fluxes  are  AU  ft  /sec,  where  A is  8/15.  In  this  report,  we  present 
numerical  results  for  three  different  U values,  2.71  ft /sec,  1.35  ft /sec, 
and  0.68  ft/sec,  which  we  designate  respectively  as  strong  flow,  medium 
flow,  and  weak  flow.  These  flow  speeds,  differing  by  factors  of  2,  have 
been  chosen  to  examine  conditions  under  which  inflow  and  outflow 
recirculation  may  arise. 
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The  portion  of  the  left  boundary,  x=0,  where  u is  positive 
requires  boundary  conditions  specifying  v,  E,  and  T.  In  conformity 
with  equations  (l6),  we  applied  the  boundary  conditions 


w = 0 , 

E = U2 /5  , 
T = 1*0 


(39a) 

(39b) 

(39c) 


for  -1*  < z < -3.  We  took  the  horizontal  domain  as  0 1 x £ 8 ft  for 

the  simulations  reported  here;  the  use  of  larger  values  of  x did  not 

significantly  change  the  results.  Our  choice  of  the  four  parameters 

L,  c.,  c , and  c for  these  numerical  simulations  was  based  on  previous 
f s w 

experience  in  the  simulation  cf  turbulent  submarine  wakes  in  a 
stratified  ocean  (Roberts  and  Piacsek,  1975;  Warn-Varnas  and  Piacsek, 
1976),  together  with  a certain  amount  of  testing.  We  adopted  the 
following  values 


T = 

±J 

0.1*  ft.  , 

(l*0a) 

cf  = 

0.5 

(U0b) 

cs 

0.1 

( l*0c ) 

c = 

0.1 

(l*0d) 

For  the  "ambient  ocean"  in  equations  (l8),  we  took  T as  60,  d as  2, 

“6  *2  2 * 

and  z very  large,  T^/z^  as  10,  Eq  as  6.7  x 10  ft  /sec  , and 

z as  2.  These  input  parameters  give  the  required  uniform  ambient 

£ 

temperature  gradient  and  an  E sufficiently  small  to  be  innocuous. 

& 

A major  difference  between  the  numerical  simulation  and  the 
proposed  experiment  is  in  the  bottom  boundary  condition  (19b)  that  the 
normal  derivative  of  u is  zero.  In  the  experiement , u vanishes  at  the 
bottom  no-slip  boundary.  This  difficulty  can  be  overcome  by  using  an 
inflow  to  the  tank  which  is  not  at  the  bottom,  and  by  modifying  NRFLC2  so 
that  the  boundary  condition  (19b)  can  optionally  be  replaced  by  zero  u. 
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The  numerical  results  for  the  three  flow  cases  are 
displayed  in  Figures  3 to  8,  which  are  photographic  reproductions 
of  the  computer  printout.  Figures  3 and  k are  for  the  case  designated 
as  strong  flow,  Figures  5 and  6 are  for  the  medium  flow  case,  and 
Figures  7 and  8 are  for  weak  flow.  The  format  is  the  same  for  each 
pair  of  figures.  In  addition.  Table  1 provides  a summary  of  significant 
numerical  values,  for  comparison  purposes. 

The  upper  panel  of  Figures  3,  5,  and  7 shows  the  mean 
temperature  solution  T(x,z),  for  the  three  inflow  speeds.  The 
contour  increment  is  5°F  in  each  case,  so  the  curves  labeled  9,  0,  1, 

2,  3,  4,  and  5 represent  Fahrenheit  temperatures  of  1*5,  50,  55 , 60, 

65,  70,  and  75.  For  strong  flow  (Figure  3)  almost  half  the  region  is 
between  1*0°F  and  1*5°F,  while  for  weak  flow  (Figure  7)  the  temperature 
distribution  is  only  slightly  disturbed  from  the  ambient  uniform 
stratification. 

The  center  panel  of  Figures  3,  5,  and  7 is  a contour  plot 

of  the  streamlines  of  the  mean  flow,  for  the  three  cases.  The  units 
2 3 

are  ft  /sec,  or  ft  /sec  per  foot  of  tank  width.  The  mean  flow  com- 
ponents u and  w are  \p  and  -\p  , where  \p  is  the  stream- function. 

^ 2 

The  contour  increment  is  0.3,  0.15,  and  0.1  ft  /sec  in  the  three  cases, 
and  the  corresponding  total  influx  and  efflux  at  the  left-hand  side 

p 

(SU/15)  is  1.1*1*,  0.72,  and  0.36  ft  /sec.  Comparison  of  the  three 
figures  indicates  that  for  strong  flow  most  of  the  cold  influx  at  the 
bottom  left  recirculates  out  at  the  top  left  (Figure  3),  while  for 
medium  flow  (Figure  5)  only  a small  part  of  the  mean  flow  recirculates. 
For  weak  flow  (Fi«rure  7)  the  stratification  stops  all  recirculation, 
and  the  single  streamline  loop  from  the  right  boundary  shows  the  effect 
of  turbulent  entrainment  by  the  bottom  Jet. 
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FIGURE  3.  Strong  Flow  Results 

The  center  panel  shows  that  most  of  the  cold  inflow 
at  the  bottom  left  recirculates  along  the  stream- 
lines and  out  at  the  top  left.  Only  one  streamline 
extends  to  the  right  boundary.  In  the  upper  panel, 
the  temperature  is  not  constant  on  streamlines 
because  of  turbulent  mixing.  The  warm  flux  from 
the  upper  right  shares  its  heat  with  the  cold 
recirculating  flux  from  the  lower  left.  Further 
results  are  shown  in  Figure  U. 
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FIGURE  1*.  Strong  Flow  Results 


These  results  supplement  those  in  Figure  3. 
The  turbulent  diffusivity  distribution  is 
obtained  from  a turbulent  kinetic  energy 
density  equation;  it  is  largest  near  the 
left  ports,  and  decreases  towards  the  right- 
hand  boundary.  The  temperature  on  the  left 
is  barely  above  Uo°F.  The  flow  through  the 
permeable  right  boundary  is  much  smaller  than 
the  imposed  flow  at  the  left,  because  of 
recirculation. 
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FIGURE  5-  Medium  Flow  Results 


In  the  center  panel,  only  one  streamline 
recirculates  from  the  bottom  left  cold  inflow 
to  the  top  left , while  three  go  out  at  the 
right  boundary.  Three  other  streamlines  bring 
warm  water  in  at  the  top  right,  and  leave  the 
region  at  the  top  left.  The  upper  panel  shows 
the  resulting  temperature  distribution;  the 
curves  0,  2,  and  U are  50°F,  60°F,  and  70°F. 
Further  results  are  shown  in  Figure  6. 
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FIGURE  5-  Medium  Flow  Results 


FIGURE  6.  Medium  Flow  Results 


These  results  supplement  those  in  Figure  5. 
The  turbulent  kinetic  energy  generated  at 
the  left  ports  is  adveeted  by  the  mean  flow, 
and  the  derived  turbulent  diffusivity  is, 
therefore,  largest  near  the  ports  and  in  the 
lower  half  of  the  region,  as  shown  in  the 
upper  panel.  In  the  center  panels,  the  hori 
zontal  fluxes  are  smaller  at  the  right  than 
at  the  left,  and  there  is  some  recirculation 
Recirculation  and  turbulent  heat  transport 
make  the  left  temperature  much  smaller  than 
the  ambient  stratification. 
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FIGURE  6.  Medium  Flow  Results 


FIGURE  7.  Weak  Flow  Results 


In  the  center  panel,  the  cold  turbulent  jet 
at  the  bottom  left  entrains  an  extra  stream- 
line and  flows  out  at  the  right;  three  stream- 
lines carry  warm  water  from  the  top  right  to 
the  top  left;  and  there  is  no  recirculation. 
The  temperature  is  almost  constant  on  the 
streamlines , and  is  little  disturbed  from  the 
ambient  stratification.  As  shown  in  the  upper 
panel,  there  is  some  turbulent  diffusion  near 
the  left  boundary  and  in  the  bottom  jet. 
Further  results  are  shown  in  Figure  8. 


UU 


••  TcMIcAAlUAt  •• 


UNIFORM  ST  RAT If  1C ATI JN.  MfcAA  OOTTOH  Jf  T AND  10*  SUCTION. 


TINT 


U.»| 


NA4INUN*  70.0  MINIMUM*  *0.0  CONTOUR  2fR0*  0.0  CONTOUR  iNCRCMfNT*  9.00  07-01-7* 

0.0  0.0  1.0  M 1.1  A.O  A.O  0.0  A.  A U 0.0 


* immiii 

a 

aa 

A 9999959999999999999999999999 

AAA 

AAA  * *444**4*4 

AA  AAAaAAAAAAA**  *• AAA  A A AAA  AAA  44 AAA A4 AAA AAA AA *4 *4 AAA 

444*4444*44444*4444444444444444*44 

9991  444A44A4 

9999  9991)999199 

9999991  1999)999)199949  99194994949949999)999999999 

9944941)  99)9949999999999999999 

999999999999999999999 

2222222222222222222222222222222222 
2222222222222  2222221222222222112 
2222222222222222  22222222222222222 

aim  mi  in  22222222222222222222 

11111111 11 111 111  mi  111 1 1 11 1111 11  mu 

immiu  11111111111111111 

uimm  iiiiiiiiiiiiu 

111111111111111 

0 0 0 000  0000  00  JOOOOOO  900000  >030 300003 3033 00j 9090000  1 1 1 

D 300000003000  000009000000000 

000000009000000 

000 0000000900 

AAAtOA  00000000000000 

*409  ****** 

******** 

******* 

**««** 

**** 

*** 


0.0  0.0  l.A  2.4  9.2  A.O  A.O  9.* 


7.2  0.0 


• • »T4c4*  t-lNtS  UUM/SI  • UNIFORM  STRAT I F I CAT  10*.  Ml  A*  MOTION  JfT  AND  13*  SUCTION.  T INF  44.  249 

NAAINviN*  0.4  44  MINIMUM*  -4.  AT  10-07  CONTOUR  UNO*  0.0  CONTOUR  INCRf Nf NT*  1. 0000-01  07-01-7* 

0.0  0.0  l.A  2. A 9.2  A.O  A.S  4.6  6.*  7.2  0.9 


luiiiaai  i mu  turn  ii  it  I u mini  in  it  u I u urn  u i tt  iiu  i ii  tin  i in  i iii  it  a ii  min  i mu  mu  in  ii 

222  4 lllilllllll 

22222 

I 222226 122  2222222222  >2222.2222222222222.22. 2222222222222222222222222222222222222222222222222222222222222222 

91  2. 2 2. 


li)44i 49 194 994)99  949  I 3j  19 1)9  144 49 4134 9)93 9 449 39  94 9999 9)99 999499 9)99)9 >99  9 4999 49999)9 9 99 


44A4NNN4**.** *****444444*4444*444*6 *****44444444444*4444 


444444444444*4444*44444 


•44444 **646*44**4444444444444 44646*44444464444444444 


4646446646*444444444444444444444444444444 


3 99  4 994  934  49  49  493  9944  33a  i 414  J3  4 93?  39)999  939  9)9  99)9)99)  99  99  9>»3))))999  9)9  9 3 49  9 9394499  9 99 


2*222  22  22222  22:222  22 


2222222222222222222222222222222222222222222222222222222222222222222222222222222 


11 111111  i 1111  iu  11 111  mini 


1111 11111111111111 111111111111 11111111 11111 1111 1111 11 11 1111111111111111  mi  11111 11111 


u 1 t 1 

U T A 

0 7 A 

u r a 

t f A 

U T A 
0 Tl 

u 0 

OTA 


UNHOAM  STRATH  1C  ATI  JN.  hCAR  93T  T JN 

4 ICO  T J£ITM  f Ilf  Ml  Mf  NT  SIMULATION. 

0 JT  T Jm  j£|  TfMlfRATURc  IS  *0  3C«4tfS 
TfM«*4ATUAf  RANOf  43  TO  03  OfiRttS 

JfT  AND  TOR 

1 FOOT  JfT 
f AHRfNMf  IT 
AT  INFINITY 

SUCTION. 

ANO  SUCTION. 

ANIO 

0.2300  * 0.  200  0 

0.  0 

0.0 

AMI  3 

0.230  0.  200 

0.0 

0.0 

Tana  1.1 

0.  0 0. 0 

3.0 

0.0 

r«  0 

T A U 
T A 0 
T A U 
T A I 
r a u 
A IU 


Mi  90.0 
*T/%  *0.0 


TT, Tl. 2TC.MI0  00.0  40.0  0.730  1000.000 

0AM0.00C*  1. 00  - 03  0.  A3 

7i«N,C*,C4,CM.*ORC  0.100  3. *33  0.103  0.103  3.900 

»**» •*TRCN»2ST RCh»| 1*33  2.03)  3.0  0.0  12  12 

T|N€.CFt.N0UTMT,DSC4lC,Ml3*t  20.0  3.00  2 4.  0 I. 

07-31-7* 


NAAN  RIANT  I MF  . OM  T(M»A4ATU«f  IS  A*. A 


3(3  99.0  47.0  90.0  71.0  09.0 

IT/S  -0.0  0.0  0.0 


FIGURE  7.  Weak  Flow  Results 


FIGURE  8.  Weak  Flow  Results 


These  results  supplement  those  in  Figure  7. 

The  turbulence  distribution  associated  with 
the  bottom  jet  is  almost  separate  from  that 
near  the  top  left  outflow  port.  The  bottom 
right  outflow  is  larger  than  the  imposed  bottom 
left  inflow,  because  of  entrainment.  Even  on 
the  left  boundary,  the  temperature  is  not  very 
different  from  the  ambient  distribution. 


The  lower  panel  o?  Figures  3 to  8 provides  a computer 
caption,  a copy  of  the  input  data,  the  computation  date,  and  the 
mean  temperature  of  the  surface  outflow  from  the  region,  calculated 
using  equation  (17).  The  three  mean  temperatures  are  1*9. 8°F,  56.1°F, 
and  66.8°F;  the  value  for  strong  flow  is  close  to  1*0°F  because  most 
of  the  surface  outflow  has  come  direct  from  the  cold  bottom  left  inflow. 

To  the  left  and  right  of  the  lower  panels  in  Figures  3,  5, 
and  7 are  plots  of  T and  u at  the  left  and  right  boundaries  of  the 
computational  domain,  together  with  the  ambient  temperature  distribution, 
all  as  functions  of  z.  These  plots  are  repeated  at  double  size  in  the 
center  panel  of  Figures  1*,  6,  and  8. 

The  ambient  temperature  increases  linearly  from  1*0°F  to 
80°F  in  each  case.  The  temperature  on  the  left  boundary  is  Uo°F  in 
the  cold  inflow  region  at  the  bottom,  and  increases  to  temperatures 
of  52.0°F,  59- 6CF,  and  69.6°F  at  the  top  left  boundary.  The  tempera- 
ture on  the  right  boundary  deviates  somewhat  from  the  ambient  tempera- 
ture, but  only  where  the  flow  is  out  of  the  computational  domain. 

The  horizontal  flow  at  the  left  boundary  is  given  by 
equation  (38),  with  the  three  maximum  flow  speeds  U differing  by 
factors  of  two.  The  horizontal  flow  at  the  right  boundary  is  more 
than  the  estimated  maximum  possible  stable  flow  jcf.  equation  ( 36 ) j , 

u = -N(z  + 2)  , (1*1) 

by  factors  of  1.1*8,  1.32,  and  1.05  in  the  three  cases.  The  estimate 
(36)  was  based  on  the  assumption  that  E is  negligibly  small  near  the 
right  boundary;  in  our  three  calculations,  this  was  not  the  case. 
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The  flow  (1*1)  corresponds  to  an  outward  flux  at  the  bottom 
right  of  0.1*1  ft  /sec,  with  an  equal  influx  at  the  top  right.  The 
actual  fluxes  at  the  right-hand  boundary  are  the  stream  function 
maxima  on  the  boundary,  respectively,  0.6l,  0.5l*,  and  0.1*3  ft  /sec. 

For  the  strong  and  medium  flow  cases , these  fluxes  are  smaller  than  the 
imposed  values  on  the  left  boundary,  and  the  remaining  flux  recirculates 
from  the  bottom  left  to  the  top  left. 

In  all  three  cases,  the  bottom  left  inflow  jets  entrain 

some  of  the  surrounding  fluid.  The  total  jet  flux,  including  this 

entrained  water,  is  the  stream-function  maximum,  respectively,  1.1*5, 

2 

O.Tl*,  and  0.1*3  ft  /sec  in  the  three  cases. 

The  turbulent  diffusivity  K(x,z)  is  displayed  in  the  upper 

panel  of  Figures  1*,  6,  and  8.  The  contour  increments  are  0.075* 

2 

0.03,  and  0.015  ft  /sec.  This  turbulent  diffusivity  distribution  is 
the  crucial  result  of  the  turbulence  modelling,  since  it  is  the  means 
by  which  the  turbulent  kinetic  energy  density  E influences  the  mean 
flow  and  temperature.  The  diffusivity  K is  largest  at  the  inflow  and 
outflow  regions  (where  it  is  proportional  to  U),  and  decreases  away 
from  the  left  boundary.  The  minimum  value  on  the  left  wall  is  signifi- 
cant, since  it  determines  the  momentum  and  heat  transport  between  the 
cold  lower  layers  moving  to  the  right  and  the  warm  upper  layers  moving 
to  the  left.  The  three  minima  are  0.250,  0.115,  and  0.0l*3  ft  /sec, 
respectively,  and  it  appears  that  for  still  weaker  flow  the  value  would 
become  so  small  that  the  two  flows  would  be  essentially  independent. 

These  numerical  results  are  summarized  in  Table  1,  for 
•.purposes  of  comparison.  The  left  boundary  fluxes  are  proportional  to 
the  flow  speed  U,  while  the  right  fluxes  are  bounded,  roughly  speaking, 
by  equation  (1*1).  The  flux  difference  recirculates  from  the  bottom 
left  inflow  to  the  top  left  outflow.  The  bottom  turbulent  jet  entrains 
some  surrounding  fluid,  which  for  the  non-recirculating  weak  flow  case 
is  drawn  in  from  infinity,  with  the  flow  near  the  surface.  The 
temperature  at  the  top  left  is  smaller  than  80°F  only  because  of 
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TABLE  1.  COMPARISON  OF  NUMERICAL 
RESULTS  FOR  THE  THREE  FLOW  CASES 


FLOW  CASE 

STRONG 

MEDIUM 

WEAK 

Maximum  Flow  Speed  U (ft/sec) 

2.68 

1.31* 

0.67 

Left  Boundary  Fluxes  (ft  /sec) 

1.1*1* 

0.72 

0.36 

Right  Boundary  Fluxes 

0.6l 

0.5l* 

0.1*3 

Recirculation  Flux 

0.83 

0.18 

0.00 

Total  Entrained  Flux 

1.1*5 

0.7l* 

0.1*3 

Temperature  at  Top  Left  (°F) 

52.00 

59.60 

69.6O 

Mean  Surface  Outflow  Temperature 

1*9.80 

56.10 

66.80 

Maximum  Diffusivity  (ft  /sec) 

0.1+67 

0.231* 

0.117 

Minimum  Diffusivity  on  Left 

— 

0.250 

0.115 

0.01*3 

I 

turbulent  heat  transport ; while  the  mean  surface  outflow  temperature 
t is  also  reduced  by  recirculation  in  the  strong  and  medium  flow  cases. 

The  diffusivity  maxima,  and  the  minima  on  the  left  boundary,  are  also 
given  in  the  table,  for  comparison  purposes. 

In  conclusion,  we  note  that  for  this  numerical  simulation 
of  our  proposed  experiment,  relatively  small  increases  in  the  flow 
speed  can  produce  recirculation,  with  a substantial  reduction  in  the 
mean  temperature  of  the  warm  surface  outflow  from  the  experimental 
tank.  The  corresponding  situation  for  an  OTP?  prototype  would  be  a 
reduction  in  the  temperature  difference  between  the  two  intakes, 
resulting  in  a plant  power  loss.  Further  research  and  experimental 
verification  is  required  to  determine  how  such  a power  loss  can  be 
avoided. 
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